options(java.parameters = "-Xmx2G")

library(r5r)
## Loading required namespace: sf
## Loading required namespace: rJava
## Please make sure you have already allocated some memory to Java by running:
##   options(java.parameters = '-Xmx2G').
## Currently, Java memory is set to -Xmx2G
library(osmextract)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
## Check the package website, https://docs.ropensci.org/osmextract/, for more details.
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.3     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(ggthemes)
library(ggspatial)
library(tigris)
## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
library(wesanderson)
library(tidytransit)

Load Networks

dir.create("networks")
## Warning in dir.create("networks"): 'networks' already exists
download.file("https://app.mecatran.com/urb/ws/feed/c2l0ZT1zbG90cmFuc2l0O2NsaWVudD1zZWxmO2V4cGlyZT07dHlwZT1ndGZzO2tleT0zZTMwMzM1OTRiMTE2NzA0N2IxNjQwNjA0ZjQwMGMzMzdiM2E1MTQ0", file.path("networks","SLOgtfs.zip"), mode = "wb", quiet=TRUE)

Street Network Data

SLO_file <- oe_match("San Luis Obispo")
## No exact match found for place = San Luis Obispo and provider = geofabrik. Best match is Manitoba. 
## Checking the other providers.
## An exact string match was found using provider = openstreetmap_fr.

Download data to networks folder and read layer of lines representing the street network

SLO_streets <- oe_read(SLO_file$url, 
                   provider = "openstreetmap_fr", 
                   download_directory = "networks", 
                   layer = "lines", 
                   quiet = TRUE) %>%
  filter(!is.na(highway))

Test Plot

ggplot(SLO_streets) +
  geom_sf()

Mask Streets to be only in the City Limits and Project onto Coordinate System

CA5_state_plane <- "+proj=lcc +lat_1=34.03333333333333 +lat_2=35.46666666666667 +lat_0=33.5 +lon_0=-118 +x_0=2000000 +y_0=500000 +ellps=GRS80 +units=m +no_defs"

SLO_city_limits <- places("California") %>%
  filter(NAME == "San Luis Obispo") %>%
  st_transform(crs = st_crs(SLO_streets))
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%
SLO_streets <- SLO_streets[SLO_city_limits,]

ggplot(SLO_streets) +
  geom_sf() +
  coord_sf(crs = CA5_state_plane)

Load School Locations and rename osm_id column to be “id” for the r5r package to read correctly

SLO_schools <- oe_read(SLO_file$url, 
                   provider = "openstreetmap_fr", 
                   download_directory = "networks", 
                   layer = "points", 
                   quiet = TRUE) %>%
  filter(str_detect(other_tags, '"amenity"=>"school"')) %>%
  st_filter(SLO_city_limits) %>%
  rename(id = osm_id)

Create grid of points (hexagons) across entire city to compare travel times to/from

grid <- st_sf(st_make_grid(SLO_city_limits, 
                           square = FALSE, 
                           n = c(100,100),
                           what = "polygons")) %>%
  st_filter(SLO_city_limits) 

colnames(grid) <- "geometry"
st_geometry(grid) <- "geometry"

grid <- grid %>%
  mutate(id = seq(1, length(grid$geometry), by=1))

grid_points <- st_centroid(grid)
## Warning in st_centroid.sf(grid): st_centroid assumes attributes are constant
## over geometries of x
ggplot(grid) +
  geom_sf() +
  geom_sf(data = grid_points, alpha = 0.5) +
  geom_sf(data = SLO_streets) +
  coord_sf(crs = CA5_state_plane) + 
  theme_map()

Set up r5r core

r5r_core <- setup_r5("networks", verbose = FALSE)
## Using cached version from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/r5r/jar/r5-v6.2-all_20210408.jar
## 
## Using cached network.dat from networks/network.dat

Calculate travel time matrix

ttm <- travel_time_matrix(r5r_core = r5r_core,
                          origins = SLO_schools,
                          destinations = grid_points,
                          mode = c("WALK", "TRANSIT"),
                          departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                          max_walk_dist = 1000,
                          max_trip_duration = 480,
                          verbose = FALSE)
## Warning in assert_points_input(destinations, "destinations"): 'destinations$id'
## forcefully cast to character.

Reformat ttm to have each row be an origin and columns for the schools

tt_wide <- ttm %>%
  pivot_wider(names_from = fromId, 
              names_prefix = "from", values_from = travel_time) %>%
  rename(id = toId) %>% 
  merge(grid) %>%
  replace(is.na(.), 999) %>%
  rowwise() %>%
  mutate(from_any = min(c_across(starts_with("from")), na.rm = TRUE))

st_geometry(tt_wide) <- "geometry"

Plot travel time to nearest school

ggplot(SLO_streets) +
  geom_sf(data = tt_wide, 
          aes(fill = from_any), 
          color = NA) +
  geom_sf(alpha = 0.5) +
  scale_fill_gradient2(low = "green", mid = "yellow", high = "red", 
                       midpoint = 30,
        name = "Transit Travel\ntime to the\nnearest school\n(minutes)") +
  coord_sf(crs = CA5_state_plane) +
  theme_map()

Create Isochrones

iso_pallete <- wes_palette("Zissou1", n = 5)

iso10min <- tt_wide[tt_wide$from_any < 11,] %>%
  st_union()

iso20min <- tt_wide[tt_wide$from_any < 21,] %>%
  st_union()

iso30min <- tt_wide[tt_wide$from_any < 31,] %>%
  st_union()

ggplot(SLO_streets) +
  geom_sf(data = iso30min, 
          aes(fill = "Area within 30 minutes"), 
          color = NA) +
  geom_sf(data = iso20min, 
          aes(fill = "Area within 20 minutes"), 
          color = NA) +
  geom_sf(data = iso10min, 
          aes(fill = "Area within 10 minutes"), 
          color = NA) +
  geom_sf(alpha = 0.5) +
  scale_fill_manual(values = c(iso_pallete[1], 
                               iso_pallete[3],
                               iso_pallete[5]),
        name = "Transit Travel\ntime to the\nnearest school\n(minutes)") +
  coord_sf(crs = CA5_state_plane) +
  theme_map()

Load transit stop locations

SLO_transit <- read_gtfs(file.path("networks", "SLOgtfs.zip"))

transit_stops <- st_as_sf(SLO_transit$stops, 
                          coords = c("stop_lon", "stop_lat"), 
                          crs =st_crs(grid))

Set up grid of points with stops data

transit_grid <- grid %>%
  mutate(num_stops = lengths(st_covers(grid, transit_stops)))

transit_points <- st_centroid(transit_grid)
## Warning in st_centroid.sf(transit_grid): st_centroid assumes attributes are
## constant over geometries of x
ggplot(transit_points) +
  geom_sf(aes(color = as.character(num_stops))) +
  scale_color_manual(values = c("gray", "cornsilk", "lightgreen"), 
                    name = "Number of\ntransit stops") +
  theme_void()

Calculate accessibility - 10 minute walk to transit stop

transit_access <- accessibility(r5r_core,
                        origins = transit_points,
                        destinations = transit_points,
                        mode = "WALK",
                        opportunities_colname = "num_stops",
                        decay_function = "step",
                        cutoffs = 11,
                        departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                        max_walk_dist = 500,
                        time_window = 60,
                        percentiles = 50,
                        verbose = FALSE) %>%
  mutate(id = as.numeric(from_id)) %>%
  merge(grid)
## Warning in assert_points_input(origins, "origins"): 'origins$id' forcefully cast
## to character.
## Warning in assert_points_input(destinations, "destinations"): 'destinations$id'
## forcefully cast to character.
st_geometry(transit_access) <- "geometry"

ggplot(transit_access) +
  geom_sf(aes(fill = accessibility), color = NA) +
  scale_fill_viridis_c(name = "Transit stops\nwithin 10-minutes\nwalk") +
  coord_sf(crs = CA5_state_plane) +
  theme_void()

Calculate access with a decay function (more continuous) with “accessibility” decreasing by one half every five minutes

transit_access2 <- accessibility(r5r_core,
                        origins = transit_points,
                        destinations = transit_points,
                        mode = "WALK",
                        opportunities_colname = "num_stops",
                        decay_function = "exponential",
                        cutoffs = 5,
                        departure_datetime = as.POSIXct("15-11-2021 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S"),
                        max_walk_dist = 500,
                        time_window = 60,
                        percentiles = 50,
                        verbose = FALSE) %>%
  mutate(id = as.numeric(from_id)) %>%
  merge(grid)
## Warning in assert_points_input(origins, "origins"): 'origins$id' forcefully cast
## to character.
## Warning in assert_points_input(destinations, "destinations"): 'destinations$id'
## forcefully cast to character.
st_geometry(transit_access2) <- "geometry"

ggplot(transit_access2) +
  geom_sf(aes(fill = accessibility), color = NA) +
  scale_fill_viridis_c(name = "Accessiblity score") +
  coord_sf(crs = CA5_state_plane) +
  theme_void()

Turn off r5r core to save memory

stop_r5(r5r_core)
## r5r_core has been successfully stopped.
rJava::.jgc(R.gc = TRUE)

Save grid of accessibility values to use in next assignment

st_write(transit_access2, 'SLO_access.geojson', append=FALSE, quiet=TRUE )
## Warning in CPL_write_ogr(obj, dsn, layer, driver,
## as.character(dataset_options), : GDAL Error 6: DeleteLayer() not supported by
## this dataset.